% Calculates statistics about price changes from the stationary 
% distribution
%
% Jon Steinsson and Emi Nakamura, July 2006+
% + additional statistics Kurtosis; interquartiles... GLB 2020
%**************************************************************

function PCStats = RigidityStatsMultiSectorCalvoPlus(Q,F,F2,alphaCalvo,...
    SectorWeights,a,rp,rad)

agridsize = size(a,1);
rpgridsize = size(rp,1);
radgridsize = size(rad,1);
Ssize = agridsize*rpgridsize*radgridsize;

NSector = size(SectorWeights,1);

rpBig = repmat(rp,[1 agridsize radgridsize]);
rpBig = permute(rpBig,[2 1 3]);
rpBig = reshape(rpBig,[Ssize 1]);

for ii = 1:NSector

    % Calculate frequency of price change
    %***************************************************************

    PriceChange1 = (rpBig ~= F(ii).F);
    PriceChange2 = (rpBig ~= F2(ii).F2);
    PCStats.freq1(ii) = sum(PriceChange1.*Q(ii).Q);
    PCStats.freq2(ii) = sum(PriceChange2.*Q(ii).Q);
    PCStats.freq(ii) = alphaCalvo(ii)*PCStats.freq1(ii) ...
        + (1-alphaCalvo(ii))*PCStats.freq2(ii);

    % Calculate fraction of price changes that are increases
    %***************************************************************

    PriceIncrease1 = (F(ii).F > rpBig);
    PriceDecrease1 = (F(ii).F < rpBig);
    PriceIncrease2 = (F2(ii).F2 > rpBig);
    PriceDecrease2 = (F2(ii).F2 < rpBig);

    PCStats.fracup1(ii) = sum(PriceIncrease1.*Q(ii).Q)/PCStats.freq1(ii);
    PCStats.fracup2(ii) = sum(PriceIncrease2.*Q(ii).Q)/PCStats.freq2(ii);
    PCStats.fracup(ii) = (alphaCalvo(ii)*PCStats.freq1(ii)/PCStats.freq(ii))*PCStats.fracup1(ii) ...
        + ((1-alphaCalvo(ii))*PCStats.freq2(ii)/PCStats.freq(ii))*PCStats.fracup2(ii);


    % Calculate absolute size of price changes, price increases and
    % price decreases
    %****************************************************************

    AbsSizePC1  = abs(F(ii).F-rpBig);
    AbsSizePC2  = abs(F2(ii).F2-rpBig);

    PCStats.sizepc1(ii) = sum(AbsSizePC1.*Q(ii).Q)/PCStats.freq1(ii);
    PCStats.sizepc2(ii) = sum(AbsSizePC2.*Q(ii).Q)/PCStats.freq2(ii);
    PCStats.sizepc(ii) = (alphaCalvo(ii)*PCStats.freq1(ii)/PCStats.freq(ii))*PCStats.sizepc1(ii) ...
        + ((1-alphaCalvo(ii))*PCStats.freq2(ii)/PCStats.freq(ii))*PCStats.sizepc2(ii);

    PCStats.sizeup1(ii) = sum(AbsSizePC1.*Q(ii).Q.*PriceIncrease1)/(PCStats.freq1(ii)*PCStats.fracup1(ii));
    PCStats.sizeup2(ii) = sum(AbsSizePC2.*Q(ii).Q.*PriceIncrease2)/(PCStats.freq2(ii)*PCStats.fracup2(ii));
    PCStats.sizeup(ii) = (alphaCalvo(ii)*PCStats.freq1(ii)/PCStats.freq(ii))*PCStats.sizeup1(ii) ...
        + ((1-alphaCalvo(ii))*PCStats.freq2(ii)/PCStats.freq(ii))*PCStats.sizeup2(ii);

    PCStats.sizedw1(ii) = sum(AbsSizePC1.*Q(ii).Q.*PriceDecrease1)/(PCStats.freq1(ii)*(1-PCStats.fracup1(ii)));
    PCStats.sizedw2(ii) = sum(AbsSizePC2.*Q(ii).Q.*PriceDecrease2)/(PCStats.freq2(ii)*(1-PCStats.fracup2(ii)));
    PCStats.sizedw(ii) = (alphaCalvo(ii)*PCStats.freq1(ii)/PCStats.freq(ii))*PCStats.sizedw1(ii) ...
        + ((1-alphaCalvo(ii))*PCStats.freq2(ii)/PCStats.freq(ii))*PCStats.sizedw2(ii);
    
    
    
    
   
weight1=alphaCalvo(ii)*Q(ii).Q;
weight2=(1-alphaCalvo(ii))*Q(ii).Q;

weight=[weight1 ; weight2];
abssize=[AbsSizePC1 ; AbsSizePC2];

weight_all=weight((abssize>.0),:)./sum(weight((abssize>.0),:));
abssize_all=abssize((abssize>.0),:);

% fprintf('verif poids=1');
% sum(weight_all);
% sum(abssize_all.*weight_all);

  
        [sortx,order] = sort(abssize_all);
        sortw = weight_all(order);

         midpoint = sum(sortw)/2;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
            if csumw(j) == midpoint
              wmed = mean(sortx([j j+1]));
            else
              wmed = sortx(j+1);
            end
  
            
        [sortx,order] = sort(abssize_all);
        sortw = weight_all(order);

         midpoint = sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
            if csumw(j) == midpoint
              w25 = mean(sortx([j j+1]));
            else
              w25= sortx(j+1);
            end
            
            [sortx,order] = sort(abssize_all);
        sortw = weight_all(order);

         midpoint = 3*sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
            if csumw(j) == midpoint
              w75 = mean(sortx([j j+1]));
            else
              w75 = sortx(j+1);
            end
  PCStats.wmed(ii)=wmed;
  PCStats.w25(ii)=w25;
  PCStats.w75(ii)=w75; 
    
          % KURTOSIS
    %****************************************************************

    SizePC1  = F(ii).F-rpBig;
    SizePC2  = F2(ii).F2-rpBig;
    m1 = sum(SizePC1.*Q(ii).Q)/PCStats.freq1(ii);
    m2 = sum(SizePC2.*Q(ii).Q)/PCStats.freq2(ii);
    
    %test=sum( ((SizePC1-m1).^4).*Q(ii).Q )
    
    test=(sum(((SizePC1-m1).^4).*Q(ii).Q)-(1-PCStats.freq1(ii))*m1^4)/PCStats.freq1(ii);
    test2=(sum(((SizePC2-m2).^4).*Q(ii).Q)-(1-PCStats.freq2(ii))*m2^4)/PCStats.freq2(ii);
    kur  = (alphaCalvo(ii)*PCStats.freq1(ii)/PCStats.freq(ii))*test ...
        + ((1-alphaCalvo(ii))*PCStats.freq2(ii)/PCStats.freq(ii))*test2;
    
    test_var=(sum(((SizePC1-m1).^2).*Q(ii).Q)-(1-PCStats.freq1(ii))*m1^2)/PCStats.freq1(ii);
    test_var2=(sum(((SizePC2-m2).^2).*Q(ii).Q)-(1-PCStats.freq2(ii))*m2^2)/PCStats.freq2(ii);
    var  = (alphaCalvo(ii)*PCStats.freq1(ii)/PCStats.freq(ii))*test_var...
        + ((1-alphaCalvo(ii))*PCStats.freq2(ii)/PCStats.freq(ii))*test_var2;
    

   % PCStats.kur(ii)=kur2;
    PCStats.var(ii)=var^2;  
 
    
    SizePC1  = F(ii).F-rpBig;
    SizePC2  = F2(ii).F2-rpBig;

    weight1=alphaCalvo(ii)*Q(ii).Q;
    weight2=(1-alphaCalvo(ii))*Q(ii).Q;

    weight=[weight1 ; weight2];
    size_all=[SizePC1 ; SizePC2];

    weight_all=weight((size_all~=.0),:)./sum(weight((size_all~=.0),:));
    size_allb=size_all((size_all~=.0),:);
    
    mean_s=sum(weight_all.*size_allb);
    var_s=sum(weight_all.*(size_allb-mean_s).^2);
    kur_test=sum(weight_all.*(size_allb-mean_s).^4)/var_s^2
    PCStats.kur(ii)=kur_test;
    
    
    % INTERQUARTILE
    
    SizePC1  = F(ii).F-rpBig;
    SizePC2  = F2(ii).F2-rpBig;

    weight1=alphaCalvo(ii)*Q(ii).Q;
    weight2=(1-alphaCalvo(ii))*Q(ii).Q;

    weight=[weight1 ; weight2];
    size_all=[SizePC1 ; SizePC2];

    weight_all=weight((size_all~=.0),:)./sum(weight((size_all~=.0),:));
    size_allb=size_all((size_all~=.0),:);

 
        [sortx,order] = sort(size_allb);
        sortw = weight_all(order);

         midpoint = sum(sortw)/2;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
            if csumw(j) == midpoint
              wmed = mean(sortx([j j+1]));
            else
              wmed = sortx(j+1);
            end
  
            
        [sortx,order] = sort(size_allb);
        sortw = weight_all(order);

         midpoint = sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
            if csumw(j) == midpoint
              w25 = mean(sortx([j j+1]));
            else
              w25= sortx(j+1);
            end
            
            [sortx,order] = sort(size_allb);
        sortw = weight_all(order);

         midpoint = 3*sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
            if csumw(j) == midpoint
              w75 = mean(sortx([j j+1]));
            else
              w75 = sortx(j+1);
            end
  
  PCStats.w25b(ii)=w25;
  PCStats.w75b(ii)=w75; 
  PCStats.interq(ii)=w75-w25; 
    
    
    
    
    
    
    
    
end

PCStats.avfreq1 = SectorWeights'*PCStats.freq1';
PCStats.avfracup1 = SectorWeights'*PCStats.fracup1';
PCStats.avsizepc1 = SectorWeights'*PCStats.sizepc1';
PCStats.avsizeup1 = SectorWeights'*PCStats.sizeup1';
PCStats.avsizedw1 = SectorWeights'*PCStats.sizedw1';

PCStats.avfreq2 = SectorWeights'*PCStats.freq2';
PCStats.avfracup2 = SectorWeights'*PCStats.fracup2';
PCStats.avsizepc2 = SectorWeights'*PCStats.sizepc2';
PCStats.avsizeup2 = SectorWeights'*PCStats.sizeup2';
PCStats.avsizedw2 = SectorWeights'*PCStats.sizedw2';

PCStats.avfreq = SectorWeights'*PCStats.freq';
PCStats.avfracup = SectorWeights'*PCStats.fracup';

PCStats.avsizepc = SectorWeights'*PCStats.sizepc';
PCStats.avsizeup = SectorWeights'*PCStats.sizeup';
PCStats.avsizedw = SectorWeights'*PCStats.sizedw';









